In [ ]:

    

import os

# Third-party
from astropy.io import ascii
from astropy import table
import astropy.coordinates as coord
import astropy.units as u
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as pl
pl.style.use('apw-notebook')
%matplotlib inline

import gary.coordinates as gc
import gary.dynamics as gd
import gary.potential as gp
from gary.observation import distance_modulus
from gary.dynamics.orbit import combine as combine_orbit
from gary.dynamics.core import combine as combine_ps

from ophiuchus import galactocentric_frame, vcirc, vlsr
import ophiuchus.potential as op
pot = op.load_potential('static_mw')


    

In [ ]:

    

catalog_data_path = "/Users/adrian/projects/globber/data/gc_catalogs/"
data_path = "/Users/adrian/projects/globber/data/ngc5897/"
plot_path = "/Users/adrian/projects/globber/figures/ngc5897"
if not os.path.exists(plot_path):
    os.makedirs(plot_path)


    

In [ ]:

    

# Global configuration stuff
cluster_name = "NGC 5897"
nsamples = 1024


    

In [ ]:

    

pm_gc_main = np.genfromtxt(os.path.join(catalog_data_path,"gl_2012_J2000.cat1.txt"), dtype=None, 
                           skip_header=2, 
                           usecols=[0,2,3,6,7,8,9,10,11,12,13],
                           names=['ngc_num','ra','dec','dist','dist_err','mu_ra','mu_ra_err',
                                  'mu_dec','mu_dec_err', 'vr', 'vr_err'])

pm_gc_main = table.Table(pm_gc_main)

go = ascii.read(os.path.join(catalog_data_path,"go97_table1.txt"))

all_gc = pm_gc_main
all_gc['name'] = np.array(["NGC {}".format(x) for x in all_gc['ngc_num']])
all_gc = table.join(all_gc, go, keys='name')

cluster = all_gc[all_gc['name'] == cluster_name]
cluster_c = coord.SkyCoord(ra=float(cluster['ra'])*u.degree,
                           dec=float(cluster['dec'])*u.degree,
                           distance=float(cluster['dist'])*u.kpc)


    

In [ ]:

    

xyz = cluster_c.transform_to(galactocentric_frame).cartesian.xyz
vxyz = gc.vhel_to_gal(cluster_c, pm=(cluster['mu_ra']*u.mas/u.yr,
                                     cluster['mu_dec']*u.mas/u.yr),
                      rv=cluster['vr']*u.km/u.s, 
                      galactocentric_frame=galactocentric_frame,
                      vcirc=vcirc, vlsr=vlsr)


    

In [ ]:

    

np.random.seed(42)
_distances = np.random.normal(cluster['dist'], cluster['dist_err'], size=nsamples)
cluster_samples_c = coord.ICRS(ra=(np.zeros(nsamples) + cluster['ra'])*u.degree,
                               dec=(np.zeros(nsamples) + cluster['dec'])*u.degree,
                               distance=_distances*u.kpc)

_mu_ras = np.random.normal(cluster['mu_ra'], cluster['mu_ra_err'], size=nsamples)
_mu_decs = np.random.normal(cluster['mu_dec'], cluster['mu_dec_err'], size=nsamples)
_vrs = np.random.normal(cluster['vr'], cluster['vr_err'], size=nsamples)

# ---
samples_xyz = cluster_samples_c.transform_to(galactocentric_frame).cartesian.xyz
samples_vxyz = gc.vhel_to_gal(cluster_samples_c, 
                              pm=(_mu_ras*u.mas/u.yr, _mu_decs*u.mas/u.yr),
                              rv=_vrs*u.km/u.s, 
                              galactocentric_frame=galactocentric_frame,
                              vcirc=vcirc, vlsr=vlsr)


    

In [ ]:

    

w0 = gd.CartesianPhaseSpacePosition(pos=xyz, vel=vxyz)
mean_orbit = pot.integrate_orbit(w0, dt=-0.5, nsteps=12000)

_w0 = gd.CartesianPhaseSpacePosition(pos=samples_xyz, vel=samples_vxyz)
all_w0 = combine_ps((w0,_w0))
orbit = pot.integrate_orbit(all_w0, dt=-0.5, nsteps=12000)


    

In [ ]:

    

pers = [orbit[:,i].pericenter().value for i in range(nsamples)] * u.kpc
apos = [orbit[:,i].apocenter().value for i in range(nsamples)] * u.kpc
eccs = (apos - pers) / (apos + pers)


    

In [ ]:

    

pers_xyz = np.zeros((3,len(pers)))
pers_xyz[0] = pers.value
pers_xyz = pers_xyz*u.kpc
mx = pot.mass_enclosed(pers_xyz)

pers_xyz = np.zeros((3,len(pers)))
pers_xyz[2] = pers.value
pers_xyz = pers_xyz*u.kpc
mz = pot.mass_enclosed(pers_xyz)

rtide_x = pers.to(u.pc) * (cluster['M'] / (3*mx))**(1/3.)
rtide_z = pers.to(u.pc) * (cluster['M'] / (3*mz))**(1/3.)
peri_rtide = np.mean(np.vstack((rtide_x, rtide_z)).value*rtide_z.unit, axis=0)

core_radius = cluster['Rc']*u.pc
R = peri_rtide / core_radius
print(R.mean())


    

In [ ]:

    

style = dict(
    color='#000000',
    histtype='step',
    linewidth=2.
)

fig,axes = pl.subplots(2,2,figsize=(5.5,5.5))

axes[0,0].hist(pers, bins=np.linspace(0,10,16), **style)
axes[0,0].set_xlabel(r"$r_p$ [kpc]")

axes[0,1].hist(apos, bins=np.linspace(5,25,16), **style)
axes[0,1].set_xlabel(r"$r_a$ [kpc]")

axes[1,0].hist(eccs, bins=np.linspace(0,1,16), **style)
axes[1,0].set_xlabel(r"$e$")

axes[1,1].hist(R, bins=np.linspace(0.5,10.,16), **style);
axes[1,1].set_xlabel(r"$\mathcal{R}$")
axes[1,1].set_yticks(np.arange(0,200+50,50))

fig.tight_layout()
fig.savefig(os.path.join(plot_path, "orbital-props.pdf"))


    

In [ ]:

    

cluster_c.galactic


    

In [ ]: